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We study the ground state phase diagram of the bilayer Heisenberg model on square lattice with 
a Bosonic RVB wave function. The wave function has the form of a Gutzwiller projected Schwinger 
Boson mean field ground state and involves two variational parameters. We find the wave function 
provides an accurate description of the system on both sides of the quantum phase transition. 
Especially, through the analysis of the spin structure factor, ground state fidelity susceptibility and 
the Binder moment ratio Q2 , a continuous transition from the antiferromagnetic ordered state to the 
quantum disordered state is found at the critical coupling of Oc = J±/J\\ ~ 2.62, in good agreement 
with the result of quantum Monte Carlo simulation. The critical exponent estimated from the finite 
size scaling analysis(l/!/ ~ 1.4) is consistent with that of the classical 3D Heisenberg universality 
class. 
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I. INTRODUCTION 



The study of quantum phase transition is a central is- 
sue in modern condensed matter physics. It is widely be- 
lieved that the Ginzburg-Landau- Wilson theory for clas- 
sical phase transition may fail to describe the quantum 
phase transition as a result of the quantum interference 
effect between classical paths(the Berry phase effect). 
Recently, the concepts of quantum order and de-confined 
quantum criticality are put forward theoretically. On the 
experimental side, the study of quantum phase transi- 
tion plays an important role in areas ranging from high- 
Tc cuprates, heavy Fermion systems, to the cold atom 
systems [!)]. 

The bilayer Heisenberg model(BHM) on square lattice 
is a standard model for the study of quantum phase tran- 
sition. With the increase of the interlayer coupling( Jl) 
over the intralayer coupling( Jy), the ground state of the 
system evolves from a state with antiferromagnetic long 
rang order to a quantum disordered state through a con- 
tinuous phase transition. Much theoretical and numerical 
efforts have been devoted to the study of this quantum 
phase transition. 

On theoretical side, perturbative calculations starting 
from both the ordered side(spin wave expansion) 0| and 
the disordered side (the bond operator expansion) [3] have 
been applied to the system. However, due to the biased 
nature of pertubative methods, none of them can give 
an accurate description of the system in the near vicin- 
ity of the quantum phase transition. The problem is also 
treated with the Schwinger Boson mean field theory 
Although the theory does predict a phase transition be- 
tween the antiferromagnetic ordered state and the quan- 
tum disordered state, the nature of the transition is in- 
correct. The mean field theory predicts a discontinues 
dimerization transition around J_l/ J\\ = 4.62 into a state 
composed of independent interlayer dimers, while in the 



real system, the intralayer correlation is nonzero for any 
finite Ji/Jy. 

On numerical side, the model is thoroughly studied 
by a variety of methods including the high temperature 
series expansion 6j and the quantum Monte Carlo simula- 
tion(Stochastic series expansion) [7, 8]. These numerical 
works confirm the existence of the quantum critical point 
around a = J±/ J ~ 2.52. The critical exponents is found 
to be consistent with that of the classical 3D Heisenberg 
universality class, indicating the irrelevance of the Berry 
phase effect in this phase transition. 

As the quantum phase transition occurs at zero tem- 
perature, it is natural to find a description of it in terms 
of an explicit ground state wave function. The variational 
approach to quantum phase transition has the virtue that 
it focus directly on the zero temperature behavior of the 
system and provides much more detailed information on 
the quantum critical behavior. In this regard, a RVB- 



type variational wave function ll|, Il2| had been applied 
to the study of the quantum phase transition in the BHM. 
The wave function is derived from Gutzwiller projection 
of Schwinger Boson mean field ground state. It is well 
known that such a RVB wave function can describe both 
the magnetic ordered and the quantum disordered state. 
Thus, it has the potential to provide an unbiased descrip- 
tion of the quantum phase transition in BHM. The same 
type of variational wave function has been successfully 
applied to the study of the single layer two-dimensional 
Heisenberg model jlll, [13] ■ However, for the BHM, the 
variational calculation in [l^ using such a wave function 
predicts a critical coupling ac = 3.51, which is a very 
bad estimate as compared to the result of numerical sim- 
ulation. A central issue to be addressed in this paper is 
to understand why the Bosonic RVB state, which works 
so well on square lattice, fails for the BHM and how to 
improve it. 

In this paper, we propose a RVB-type variational 
wave function with two variational parameters for the 
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BHM. Similar to IJ], our wave function is derived from 
Gutzwiller projected Schwinger Boson mean field state. 
However, in our theory the intralayer RVB pairing and 
interlayer RVB pairing are treated as two independent 
variational parameters, rather than been determined by 
mean field self-consistent equations. We find our vari- 
ational wave function provides an accurate description 
of the quantum phase transition of the BHM. We find 
the the transition is continuous. By analyzing the spin 
structure factor, ground state fidelity susceptibility and 
Binder moment ratio Q2, the critical coupling strength 
is estimated to be ac ~ 2.62, in good agreement with 
those determined from the numerical simulation. The 
critical exponent estimated from the scaling analysis of 
the Q2 data is also consistent wit that of the classical 3D 
Hciscnberg universality class. Our result indicates that 
the Bosonic RVB wave function derived from Gutzwiller 
projection of the Schwinger Boson mean field state can 
provides accurate description of the quantum phase tran- 
sition in quantum antiferromagnets. We also find that 
the failure of the Schwinger Boson mean field theory orig- 
inates from the overestimation of the tendency to form 
interlayer dimers, which is again caused by the relaxation 
of the no double occupancy constraint. 

The paper is organized as follows. In section H, we 
introduce the BHM and the Bosonic RVB wave function. 
In section HI, we present the numerical method to do 
calculation on such wave functions. In section IV, we 
present the numerical results and determine the critical 
point of the phase transition by analyzing the results of fi- 
delity susceptibility and Binder moment ratio. In section 
V, we present a discussion on related issues and conclude 
this paper. 



II. THE BILAYER HEISENBERG MODEL AND 
THE RVB-TYPE VARIATIONAL WAVE 
FUNCTION 

The model(BHM) we study in this paper is given by 

H = j|| sr-s; + j^^si.s?, (1) 

where Sf denotes the spin operator at site i of layer 
/x(= 1,2). X]<ij> uieans the summation over nearest- 
neighboring sites on the square lattice of each layer. 
a = J±/J\\ is the only dimensionless parameter of the 
model. When a = 0, the model describes two decoupled 
two-dimensional Heisenberg model, each of which are 
antiferromagnetic ordered at zero temperature. When 
a — >■ cxD, the system reduces to N decoupled interlayer 
dimers and the system is in a trivial quantum disordered 
state. A continuous quantum phase transition connects 
these two limits. Earlier numerical simulation shows that 
the phase transition occurs around ac = 2.520, Q. 



The Bosonic RVB wave function we will adopt in this 
study is made of coherent superposition of spin singlet 
configurations on the lattice and can be written as 

N/2 

|RVB)= ^ A{{tk,jk})Y[S{tk,jk), (2) 

{ik,jk} k=l 

in which S{ik,jk) = -^{tikijk - U^tjJ denotes the 
spin singlet pair between site and jk- ^({*fc,Jfc}) are 
the coefficients of the coherent superposition. In our 
case, A{{ik, jk}) can be written in a factorizeable form 

M{ik,jk}) = nfi?a^fcjk 

The wave function Eq.(2) can be used directly as vari- 
ational state for quantum antiferromagnet. A more ef- 
ficient and intuitively more attractive way to generate 
the RVB wave function is by Gutzwiller projection of 
Schwinger Boson mean field state. This approach is used 
to study two-dimensional Heisenberg model and is proved 
to be very successful. However, direct application of the 
approach to the BHM results in unsatisfactory results. 

Here, we will adopt the form of the Gutzwiller pro- 
jected Schwinger Boson mean field state, but regard 
the mean field order parameters (intralayer and interlayer 
RVB pairing amplitudes) as free variational parameters, 
rather than been determined from the mean field self- 
consistent equations. The reason for such a choice is as 
follows. In the mean field treatment, the no double occu- 
pancy constraint is relaxed. As a result, the quantitative 
prediction of the mean field theory is not reliable. For 
example, the mean field equation predicts an un-physical 
dimerization transition for BHM at a ~ 4.62, whose ori- 
gin can be traced back to the overestimation the tendency 
to form interlayer dimers, which is again related to the 
relaxation of the local constraint. 

In the Schwinger Boson representationlij, the spin op- 
erator is written as 



(3) 



a, 13=1,2 



in which ba is a Boson operator, a is the Pauli matrix. 
Eq.(3) is a faithful representation of the spin algebra pro- 
vided that the Bosonic particle satisfy the no double oc- 
cupancy constraint 



(4) 



The BHM written in terms of the Schwinger Boson op- 
erators reads 



H = - j|| J2 

<i,j>.,fj. i 
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in which 
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bis..lbi,2,t) 



(6) 



denote the intralayer and interlayer RVB pairing opera- 
tor, Ui^fj, = Y^a=t,i^\t^,abi,p.,a- The Largrange multiplier 
Ai,^ is introduced to keep track of the local constraint. 

In the mean field theory, we treat A^,^ = A as a con- 
stant and decouple the interaction term using the follow- 
ing mean field order parameters A|| = (A- ^) = i^i.j) 
and Aj^ = (Ai). The mean field Hamiltonian reads(up 
to a constant) 



Hmf = -J||A,| (^'^ 



-J^A^5](Al + A„ 



A^O 



Ay^rij^^. 



(7) 



The mean field ground state of Eq.(7) reads 
\G) = exp[ ^ a.^,,. - ^l.^^^l.t)] o)(8) 



in which |0) denotes the vacuum of the Schwinger Bo- 
son, a^^jv represents the RVB amplitude between site 
i in fi layer and site j in v layer. As a result of the 
bipartite nature of the system, the RVB amplitude is 
nonzero only for sites belonging to different sublattices. 
Thus for ^ = V, ai^ji, is nonzero only when i,j have 
different parity, while for fi ^ v the reverse is true. 
The intralayer and interlayer RVB amplitudes are given 
by(aiiji = ai2j2,ajij2 = aj2ji by symmetry) 



ail,7"l 



aii,72 



-E 

k 



[S_{k) + ri{k)] exp(ifc • 



?7(fc)] exp(i/c • r^j). 



(9) 



in which 



m 



ci7(fc) 



C2 



1 + v/1- (ci7(fc) + C2)2 

Ci7(fc) - C2 



1 + Vl-(ci7(fc)-C2)2' 

JiAi/V2A, 7(fc) 



here ci = 4J|| A||/\/2A, C2 
(cos(fc2:) -I- cos(fcj^))/2. 

The Bosonic RVB wave function adopted in this study 
is given by Gutzwiller projection of the mean field ground 
state into the physical subspace satisfying the local con- 
straint. 



\G)=Pg[ 



''itj.,i"jv,t 



N/2 



(10) 



Here Pq denotes the Gutzwiller projection and N is the 
number of lattice sites. The mean field ground state con- 
tains two dimensionless parameters, namely Ci and C2. In 
the mean field theory, both of them are determined by 
the mean field self-consistent equations. Here we regard 
them as two independent variational parameters. This is 
the key difference between our theory and that of [l3| • 

The proposed wave function Eq.(lO) can describe both 
the magnetic ordered and the quantum disordered state. 
As can be seen from Eq.(9), as ci -t- C2 — ^ 1, both anji 
and flii j-2 becomes long ranged and the wave function de- 
scribes a state with antiferromagnetic long range order. 
On the other hand, when C1+C2 deviates from 1, the RVB 
amplitudes an ji and a^i j2 become short ranged and the 
corresponding wave function describes a quantum disor- 
dered state. In fact, ci + C2 = 1 is nothing but the Bose 
condensation condition in the mean field theory. 

On general grounds, we expect the interlayer pairing 
C2 to increase with a and the intralayer pairing ci to 
decrease with a. The transition between the antiferro- 
magnetic ordered state and the quantum disordered state 
is signaled by the deviation of ci -I- C2 from 1. These ex- 
pectations are confirmed in the numerical calculation. 



III. THE NUMERICAL TECHNIQUES 

The Bosonic RVB wave function Eq.(lO) can be stud- 
ied by the standard loop gas Monte Carlo algorithm [l^. 
In this algorithm, the calculation of expectation value of 
a physical quantity ^(for example the energy) is done as 
follows 



{G\A\G) _ ^;v'7'(7iy 

(G|G) 



h i Ah') 
(7I7') 



(11) 



Here I7) denotes the valence bond basis vector and is 
given by I7) = Il^ij),^^ S{i, j). tp^ is the corresponding 



amplitude and is given by ■(/'t, = H 



(ij)e7 "-^J 



The overlap 



between two valence bond basis vectors I7) and I7') can 
be graphically interpreted as a loop gas on the lattice by 
fusing the valence bonds in the two basis vectors. It is 
easy to show that (7I7') = 2^^ , where Nj^ is the number 
of loops in the transition graph between I7) and I7'). 

As the system is bipartite, the RVB amplitude anji 
and aii,j2 are in fact positive definite and the wave func- 
tion Eq.(lO) satisfy the Marshall sign rule[l3|- For this 



as 



reason, we can interpret 1^(7,7') — : — ; — , , 

z^^,-,/ '/-;'/'-,' (7i7'> 

a normalized probability in the space of loop gas and 
can draw samples on it with the standard Monte Carlo 



method. The calculation of 
and the result reads 



(7l^l7') 
(7I7') 



is easy for A = Si- S, 



(Si • Sj 



|, j,j G same loop, different sublattices; 
i,j e same loop, same sublattice; (12) 
otherwise. 
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Thus both the energy and spin structure factor can be 
easily calculated with the standard Monte Carlo proce- 
dure in the loop gas space. 

To determine the optimal value of the variational pa- 
rameter ci and C2, we calculate the expectation value of 
the energy and of its gradients in the parameter space 
(ci,C2) on a finite lattice. It is useful to note that the 
gradients of energy can be directly simulated by the loop 
gas Monte Carlo method also. Its expression is given by 

dE{ci,C2) ^ / ( 91nfl,,A {l\H\y) \ 

\(»j)e7.7' ' ' i 

where () l denotes average over the loop gas configura- 
tions with the weight 14^(7,7'). 

We have used 10^ samples to calculate the energy and 
its gradients to determine the optimized values for ci and 
C2 . The boundary condition of the finite lattice is set to 
be periodic in both directions. The calculation is done 
on a lattice with size up to 20 x 20 x 2, at which wc find 
the critical coupling converges to ac ~ 2.62. 



IV. THE NUMERICAL RESULTS 

The optimized value for the parameter ci and C2 as 
functions of the coupling constant a are shown in Fig.l. 
As a increases, the interlayer RVB pairing strength C2 
increases at the expense of the intralayer RVB pairing 
strength ci. The result is obtained on a 20 x 20 x 2 
lattice. It is found that the optimized values deviate sig- 
nificantly from the mean field predictions, especially for 
large value of a. For example, the mean field theory pre- 
dicts that A_L would reach twice the value of Ay around 
a = 4. However, the variational theory predicts that Aj^ 
is slightly smaller than Ay around a = 4. Thus, the 
mean field theory overestimates greatly the tendency to 
form interlayer dimer at large a. 

To better understand the evolution of the variational 
parameters as functions of a, we plot the value of the 
a = ci -I- C2 as a function of a in Fig. 2. As we have 
shown above, the quantum phase transition between the 
magnetic ordered state and the quantum disordered state 
in our variational theory is solely controlled by the value 
of a. The value of a is seen to deviate from unity around 
2.6, at which the Bose condensate of the spinon is gone. 

To further characterize the quantum phase transition 
and determine the value of the critical coupling ac, we 
study the following three kinds of quantities: the spin 
structure factor at the ordering wave vector, the fidelity 
susceptibility of the ground state and the Binder moment 
ratio Q2- 




0.0 — I — ' — I — ' — I — ' — I — ' — I — ' — I — ' — I — ' — I — ' — r- 
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FIG. 1: The optimized value for the intralayer and interlayer 
RVB pairing strength ci and C2 as functions of the coupling 
constant a. 
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FIG. 2: The optimized value for the variational parameter 
a — ci + C2, which controls the order-disorder transition of 
BHM. 

A. Spin Structure Factor 

For a finite system, the spontaneous magnetization can 
be defined in a spin rotational invariant way as the square 
root of the spin structure factor at the magnetic Bragg 
vector. For BHM, the Bragg vector is Q = (vr, tt, tt). The 
spin structure factor is defined as 

^(9) = ^ E(^^ ■ e^P(*9 ■ ^"'^^O (14) 

For q ~ Q, we have 

= NSiQ) = ^(-1)'-^"(S. • S,) (15) 

i,3 
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In the quantum disordered state, as the spin correlation 
length is finite, S{Q) is of order one. However, in the 
magnetic ordered state, S{Q) should scale like N and 
thus M is an extensive quantity. 

The result of the spin structure factor for a 20 x 20 x 2 
system is shown in Fig. 3. An order-disorder transition 
can be seen around 2.5. However, the signature of phase 
transition in the spin structure factor is not sharp enough 
for an accurate determination of the critical coupling 
strength. The transition is rounded into a crossover as a 
result of the finite size effect. For this reason, we need 
some other quantities that are more sensitive to the tran- 
sition to determine the critical coupling. 



the overlap between variational ground states for nearby 
values of a. The overlap between the Bosonic RVB states 
is calculated in the following way. 

(*l*> E^,yV';V'7'(7i7') ■ ^ ' 

In our calculation, we have set 6a — 0.01. The result 
for the fidelity susceptibility for systems of several sizes 
are shown in Fig. 4. A pronounced peak appears around 
a = 2.6. Fig. 5 shows the peak position extracted from 
Fig. 4 as a function of the system size L. It is found that 
the peak position converges rapidly to its thermodynamic 
limit value ac ~ 2.62 when L > 10. 




a 



FIG. 3: The spin structure factor at the antiferromagnetic 
ordering wave vector as a function of a for a system with 
L = 20. 



FIG. 4: The fidelity susceptibility for system of size L = 
6, 8, 10, 14, 16, 18 and 20 as functions of a. The finite differ- 
ence used to calculate the fidelity susceptibility is 5a — 0.01. 



B. Fidelity Susceptibility 

The concept of fidelity susceptibility is introduced to 
describe the sensitivity of the ground state to the varia- 
tion of the parameters in HamiltonianQ and is expected 
to reach its maximum at the critical coupling of a quan- 
tum phase transition, where the ground state is the most 
susceptible to the variation of the controlling parameters 
of the phase transition. The fidelity susceptibility is de- 
fined in the following manner for a system with only one 
parameter a, 



Yf = —2 lim 



ln\0{a,6a)\ 
{5af ' 



(16) 



in which 0{a,6a) — {'i>a\'i'a+Sa) denotes the overlap be- 
tween the normalized ground state vector for parameter 
value a and a + 6a. 

In our variational theory, the fidelity susceptibility can 
be calculated directly. We first fit the optimized varia- 
tional parameters as functions of a and then calculate 



C. Binder Moment Ratio Q2 



To confirm the result derived from the fidelity suscep- 
tibility, we calculate the Binder moment ratio 
The Binder moment ratio Q2 is a dimensionless quantity 
defined in the following manner. 



Q2 



(18) 



in which Sq = J2i j -S^i ' Sj. Note our definition of 

Q2 is slightly different from the standard one in that it 
is defined in a spin rotational invariant way, while in the 
standard definition only the z-component of the moment 
is used. The Binder moment ratio is very useful in the 
analysis of the critical properties as it is universal near 
the critical point. More specifically, it can be expressed 
as a universal scaling function of tL^/'^ , where t = (a— etc) 
and u is the critical exponent for correlation length. 



FIG. 5: The critical coupling strength olc. determined from 
the peak position of the fidelity susceptibility. The peak po- 
sition is seen to converge rapidly to its thermodynamic limit 
when L > 10. The error bars are smaller than the size of the 
symbols. 
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FIG. 6: The Binder moment ratio Q2 for L 
20 as functions of a. 



14, 16, 18 and 



The results of Q2 for system with L = 14, 16, 18 and 20 
are shown in Fig. 6. It is found that all curves cross with 
each other at approximately the same value of a, in accor- 
dance with the scaling hypothesis. The estimated value 
of the critical coupling strength is 2.62, in good agrement 
with that estimated from the fidelity susceptibility data. 
The Qi value at the crossing point is found to be approx- 
imately 1.23, close but smaller than the value(1.29) es- 
timated from the quantum Monte Carlo simulation with 
the standard definition of Q2- Such a difference may be 
caused by the difference in the definitions of Q^- 

Fig. 7 shows the scaling of the Qi data with the scaling 
form Q2 = M{tL^^'^), where t = a — ac and v is the ex- 



FIG. 7: Scaling of the Binder moment ratio Q2 data for L = 
14, 16, 18 and 20. Here t ^ a - Oc 



ponent for the correlation length. The best fit is reached 
by ac ~ 2.62 and l/i/ « 1.4. The critical exponent so 
obtained ly ~ 0.714 is thus quite close to the result of the 
quantum Monte Carlo simulation. 



V. CONCLUSION 

In this work, we proposed a Bosonic RVB wave func- 
tion with the form of the Gutzwiller projected Schwinger 
Boson mean field ground state for the BHM. We find the 
proposed wave function predicts a continuous phase tran- 
sition between the antiferromagnetic ordered state and 
the quantum disordered state. To determine the critical 
coupling strength, we have calculated the spin structure 
factor, the fidelity susceptibility and the Binder moment 
ratio Q2- Through finite size scaling analysis of the lat- 
ter two quantities, we find the critical coupling to be 
given by ac ~ 2.62, in good agreement with the quan- 
tum Monte Carlo simulation results. The scaling analysis 
of Q2 also provides an estimate of the correlation length 
critical exponent(l/j/ « 1.4), which is also in good agree- 
ment with the result of quantum Monte Carlo simulation. 
We find the intralayer correlation is quite large at the 
phase transition point and it dominates over the inter- 
layer correlation for a twice as large the critical coupling 
strength. Thus, the phase transition has nothing to do 
with the dimerization instability. 

Our work indicates that the Bosonic RVB wave func- 
tion derived from Gutzwiller projection of the Schwinger 
Boson mean field ground state has the potential to cap- 
ture the physics of quantum phase transition with high 
accuracy. The failure of it in previous variational study 
14| can be attributed to the weakness of the mean field 
theory, which overestimate the tendency of the system to 
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form intcrlaycr dimcrs. Such a ovcrcstimation is closely 
related to the relaxation of the local constraint in the 
mean field treatment, which prohibit multiple occupation 
of dimer on a given bond, even if the mean field theory 
points to the tendency of Bose condensation of of such 
interlayer dimers. The same instability also cause the 
failure of the mean field theory itself for large a. Hence, 
the form the ground state predicted by the mean field the- 
ory is correct, however, the quantitative relation between 
mean field order parameters is less meaningful. The local 
constraint is thus indispensable for a correct description 
of the quantum antiferromagnet with the Bosonic RVB 
state. 

In this work, we have proved the usefulness of the varia- 
tional approach to the quantum phase transition in BHM. 
However, a more detailed study of the critical behav- 
ior and the excitation spectrum around the critical point 
is obviously needed to further characterize the quantum 
critical point in this system. We will leave this task to 
future investigations. 
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